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APPENDIX 



/*******************^ 

void sobol(int n) 

This is the function sobol for generating Sobol points. It returns the 
n-th d-dimensional Sobol point. (The point (0,0, 0) is skipped.) 

The point is implicitly returned through the array x. The function 
sobseq from Numerical Recipes, 1992 was used as a basis, but 
changed significantly to accommodate the parallel distributed approach. 
The function sobol can generate Sobol points skipping an initial 
part of the sequence. The constant MAXDIM, see Numerical Recipes, 
is extended to 360 which required adding more initializing data to 
the arrays ip (the primitive polynomials) , 

mdeg (their degrees) , and iv (the initial direction numbers) . The 
polynomial x is used to generate the first coordinate of Sobol points. 

***********************************^ 

B 

.llinclude "nrutil.h" 
Sfdefine MAXBIT 30 
Mfdefine MAXDIM 3 60 

i y 

^extern int d; /* actual dimension of the points */ 

Jxtern double *x; /* This returns implicitly the n-th Sobol point in x */ 

UJoid sobol (int n) 

,{ 

p int j ,k, 1; 

unsigned long i,im,ipp; 
j jj static double fac; 

U static unsigned long in, ix [MAXDIM+1] , *iu [MAXBIT+1] ; 
%D static unsigned long mdeg [MAXDIM+1] ={ 0 , MAXBIT, 1 , 2 , 3 , 3 , 4 , 4 , 

5,5,5,5,5,5, 

m 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7, 
/*16*/ 8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8, 

9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9, 
9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9, 
9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9, 
/*2 0*/ 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 

10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 
10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 

11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 
11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11, 
11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11, 
11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 
11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11, 
11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 
11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 
11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 
11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 

12, 12, 12, 12, 12, 12, 12,12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 
12,12,12 

}; 
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/* The first coordinate is set to MAXBIT, but it is never used */ 

static unsigned long ip [MAXDIM+1] = { 0 , MAXBIT, 0,1,1,2,1,4, 
/*deg 5*/ 2,13,7,14,11,4, 
/*deg 6*/ 1,16,13,22,19,25, 

/*deg 7*/ 1,32,4, 8,7,56,14,28,19,50,21,42,31, 62,37,41,55,59, 

/*deg 8*/ 14, 56, 21, 22, 3,8, 47, 49, 50, 52, 67 , 70 , 84 , 97 , 103 , 115 , 122 , 

/*deg 9*/ 8,13 ,16 ,22 ,25,44 ,47 ,52 ,55 ,59 ,62 ,67 ,74 ,81 ,82 ,87 ,91 ,94, 
103, 104, 109, 122, 124, 13 7, 13 8, 143, 145, 152, 157, 167, 173, 176, 181, 182, 185, 
191, 194, 199,218,22 0,2 2 7,22 9,23 0,234,236,241,244, 253, 
/*deg 10*/ 4, 13, 19,22, 50, 55, 64, 69, 98, 107, 115, 121, 127, 134, 140, 145, 152, 

158, 161, 171, 181, 194, 199,2 03,2 08,22 7,242,2 51,253,265,2 66,2 74,2 83,2 89,2 95, 
3 01,316,319,324,346,352,361,367,382,395,3 98,400,412,419,422,426,428,433, 
446,4 54,457,4 72,4 93, 505, 508, 
/*degll*/ 2,11,21,22,35,49,50,56,61,70,74,79,84,88,103,104,112,115,117,122,134, 
13 7, 146, 14 8, 157, 158, 162, 164, 168, 173, 185, 186, 191, 193, 199,213,214, 22 0, 22 7, 
236,242,251,256,259,265,266,276,292,304,310,316,319,322,328, 334,339,341, 
345, 346,3 62,3 67,3 72,3 75,3 76,3 81,3 85,3 88,3 92,4 0 9,415,416,421,42 8,431,434, 
43 9,446,451,453,457,458,471,475,478,484,493,4 94,4 99,502,517,518,524,52 7, 
555, 560, 565, 569, 578, 5 80, 5 87, 5 89, 590, 601, 607, 611, 614, 617, 618, 62 5, 62 8, 635, 
641,64 7,654,659,662,672,675,682,684,689,695,696,713,719,724,73 3,734,74 0, 
^ 747, 74 9, 752 , 755, 762 , 770, 782 , 784 , 787 , 789 , 793 , 7 96 , 803 , 805 , 810 , 815 , 824 ,829, 

U 83 0, 832 , 841, 84 7, 849, 861, 871, 878, 889, 892, 901, 90 8, 92 0, 923, 942, 949, 950, 954 , 

yj . 961, 968, 971, 973, 979, 982, 986, 998, 1001, 1010, 1012, 

f|*degl2*/ 41 , 52 , 61 , 62 , 76 , 104 , 117 , 131 , 143 , 145 , 157 , 167 , 171, 176 , 181, 194 , 217 , 236 , 
^ 2 3 9,262,283,286,307}; 

D static unsigned long iv [MAXDIM*MAXBIT+1] = { 
Si 0, /* 4 0 elements in a row */ 

[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, 
^,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, 
"f, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1,1, 1,1, 1,1, 
^1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, 

y|, 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, 
ill, 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, 

r^l, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1,1, 
%, 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, 

! % 1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1,1, 1,1, 1,1, 1,1,1, 1,1, 1,1,1, /*end 1 

W, 3, 1,3, 1,3, 1,3, 3,1, 3, 1,3, 1,3, 1,1, 3, 1,3, 1,3, 1,3,3, 1,3, 1,3, 1,3, 1,1,3, 1,3, 1,3, 1,3, 
1,3,1,3,1,3,1,3,3,1,3,1,3,1,3,1,1,3,1,3,1,3,1,3,3,1,3,1,3,1,3,1,1,3,1,3,1,3,1,3, 
1,3,1,3,1,3,1,3,3,1,3,1,3,1,3,1,1,3,1,3,1,3,1,3,3,1,3,1,3,1,3,1,1,3,1,3,1,3,1,3, 
1,3,1,3,1,3,1,3,3,1,3,1,3,1,3,1,1,3,1,3,1,3,1,3,3,1,3,1,3,1,3,1,1,3,1,3,1,3,1,3, 
1,3,1,3,1,3,1,3,3,1,3,1,3,1,3,1,1,3,1,3,1,3,1,3,3,1,3,1,3,1,3,1,1,3,1,3,1,3,1,3, 
1,3,1,3,1,3,1,3,3,1,3,1,3,1,3,1,1,3,1,3,1,3,1,3,3,1,3,1,3,1,3,1,1,3,1,3,1,3,1,3, 
1,3,1,3,1,3,1,3,3,1,3,1,3,1,3,1,1,3,1,3,1,3,1,3,3,1,3,1,3,1,3,1,1,3,1,3,1,3,1,3, 
1,3,1,3,1,3,1,3,3,1,3,1,3,1,3,1,1,3,1,3,1,3,1,3,3,1,3,1,3,1,3,1,1,3,1,3,1,3,1,3, . 
1,3, 1,3, 1,3, 1,3, 3, 1,3, 1,3, 1,3, 1,1, 3, 1,3, 1,3, 1,3, 3, 1,3, 1,3, 1,3,1, 1,3, 1,3, 1,3, 1,3, /*end 2 

1,5,7,7,5,1,3,3,7,5,5,7,7,1,3,3,7,5,1,1,5,3,3,1,7,5,1,3,3,7,5,1,1,5,7,7,5,1,3,3, 

1,5,7,7,5,1,3,3,7,5,5,7,7,1,3,3,7,5,1,1,5,3,3,1,7,5,1,3,3,7,5,1,1,5,7,7,5,1,3,3, 

1,5,7,7,5,1,3,3,7,5,5,7,7,1,3,3,7,5,1,1,5,3,3,1,7,5,1,3,3,7,5,1,1,5,7,7,5,1,3,3, 

1,5,7,7,5,1,3,3,7,5,5,7,7,1,3,3,7,5,1,1,5,3,3,1,7,5,1,3,3,7,5,1,1,5,7,7,5,1,3,3, 

1,5,7,7,5,1,3,3,7,5,5,7,7,1,3,3,7,5,1,1,5,3,3,1,7,5,1,3,3,7,5,1,1,5,7,7,5,1,3,3, 

1,5,7,7,5,1,3,3,7,5,5,7,7,1,3,3,7,5,1,1,5,3,3,1,7,5,1,3,3,7,5,1,1,5,7,7,5,1,3,3, 

1,5,7,7,5,1,3,3,7,5,5,7,7,1,3,3,7,5,1,1,5,3,3,1,7,5,1,3,3,7,5,1,1,5,7,7,5,1,3,3, 

1,5,7,7,5,1,3,3,7,5,5,7,7,1,3,3,7,5,1,1,5,3,3,1,7,5,1,3,3,7,5,1,1,5,7,7,5,1,3,3, 

1,5,7,7,5,1,3,3,7,5,5,7,7,1,3,3,7,5,1,1,5,3,3,1,7,5,1,3,3,7,5,1,1,5,7,7,5,1,3,3, /*end 3 

*/ 

1,15,11,5,3,1,7,9,13,11,1,3,7,9,5,13,13,11,3,15,5,3,15,7,9,13,9,1,11,7,5,15,1,15,11,5,3,1,7, 
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1,15,11,5,3,1,7,9,13,11,1,3,7,9,5,13,13,11,3,15,5,3,15,7,9,13,9,1,11,7,5,15,1,15,11,5,3,1,7,9, 
1,15,11,5,3,1,7,9,13,11,1,3,7,9,5,13,13,11,3,15,5,3,15,7,9,13,9,1,11,7,5,15,1,15,11,5,3,1,7,9, 
1,15,11,5,3,1,7,9,13,11,1,3,7,9,5,13,13,11,3,15,5,3,15,7,9,13,9,1,11,7,5,15,1,15,11,5,3,1,7,9, 
1,15,11,5,3,1,7,9,13,11,1,3,7,9,5,13,13,11,3,15,5,3,15,7,9,13,9,1,11,7,5,15,1,15,11,5,3,1,7,9, 
1,15,11,5,3,1,7,9,13,11,1,3,7,9,5,13,13,11,3,15,5,3,15,7,9,13,9,1,11,7,5,15,1,15,11,5,3,1,7,9, 
1,15,11,5,3,1,7,9,13,11,1,3,7,9,5,13,13,11,3,15,5,3,15,7,9,13,9,1,11,7,5,15,1,15,11,5,3,1,7,9, 
1,15,11,5,3,1,7,9,13,11,1,3,7,9,5,13,13,11,3,15,5,3,15,7,9,13,9,1,11,7,5,15,1,15, 11,5,3,1,7,9, 
1,15,11,5,3, 1,7,9,13,11,1,3,7,9,5,13,13,11,3,15,5,3,15,7,9,13,9,1,11,7,5,15,1,15, 11,5,3,1,7,9, 



1,17,13,7,15,9,31,9,3,27,15,29,21,23,19,11,25,7,13,17,1,25,29,3,31,11,5,23,27,19,21,5,1,17,13, 
7,15,9,31,9, 

1,17,13,7,15,9,31,9,3,27,15,29,21,23,19,11,25,7,13,17,1,25,29,3,31,11,5,23,27,19,21,5,1,17,13, 
7,15,9,31,9, 

1,17,13,7,15,9,31,9,3,27,15,29,21,23,19,11,25,7,13,17,1,25,29,3,31,11,5,23,27,19,21,5,1,17,13, 
7,15,9,31,9, 

1,17,13,7,15,9,31,9,3,27,15,29,21,23,19,11,25,7,13,17,1,25,29,3,31,11,5,23,27,19,21,5,1,17,13, 
7,15,9,31,9, 

1,17,13,7,15,9,31,9,3,27,15,29,21,23,19,11,25,7,13,17,1,25,29,3,31,11,5,23,27,19,21,5,1,17,13, 
7,15,9,31,9, 

1,17,13,7,15,9,31,9,3,27,15,29,21,23,19,11,25,7,13,17,1,25,29,3,31,11,5,23,27,19,21,5,1,17,13, 
J, 15, 9, 31, 9, 

\M, 17, 13, 7, 15, 9, 31, 9, 3, 27, 15, 29, 2 1,23, 19, 11, 25, 7, 13, 17, 1,25 ,29, 3, 31, 11, 5 ,23, 27, 19, 21, 5 ,1,17, 13, 
ig, 15, 9, 31, 9, 

rt, 17, 13, 7, 15, 9, 31, 9,3,27,15,29,21,23,19,11,25,7,13,17,1,25,29,3,31,11,5,23,27,19,21,5,1,17,13, 
£? u 15, 9, 31, 9, 

1^,17,13,7,15,9,31,9,3,27,15,29,21,23,19,11,25,7,13,17,1,25,29,3,31,11,5,23,27,19,21,5,1,17,13, 
S, 15, 9, 31, 9, /*end 5 */ 

1^,51,61,43,51,59,47,57,35,53,19,51,61,37,33,7,5,11,39,63,27,17,15,23,29,3,21,13,31,25,9,49,33, 
f£9, 29,11, 19,27,15,25, 

% 51, 61, 43, 51, 59, 47, 57, 35, 53, 19, 51, 61, 37, 33, 7, 5, 11, 39, 63, 27, 17, 15, 23, 29, 3, 2 1,13, 31, 25, 9, 4 9, 33, 
^19, 29, 11, 19, 27, 15, 25, 

S, 51, 61, 43, 5 1,59, 47, 57, 3 5, 53, 19, 51, 61, 37 ,33, 7, 5, 11, 39, 63, 27, 17,15,23,29,3,21,13,31,25,9,49,33, 
f|J9,29, 11, 19, 27, 15, 25, 

k£, 51,61,43,51,59,47,57,35,53,19,51,61,37,33,7,5,11,39,63,27, 17,15,23,29,3,21,13,31,25,9,49,33, 
1:9, 29, 11, 19,27, 15,25, 

% 51, 61, 43, 51, 59, 47, 57, 35, 53, 19, 51, 61, 37, 33, 7, 5, 11 ,39, 63, 27, 17,15,23,29,3,21,13,31,25,9,49,33, 
y£9,29, 11, 19,27, 15,25, 

m, 51,61,43,51,59,47, 57,35,53,19,51,61,37,33,7,5,11,39,63,27, 17,15,23,29,3,21,13,31,25,9,49,33, 
"19,29,11,19,27,15,25, 

1,51,61,43,51,59,47,57,3 5,53,19,51,61,37,33,7,5,11,39,63,27, 17,15,23,29,3,21,13,31,25,9,4 9,33, 
19,29,11,19,27,15,25, 

1,51,61,43,51,59,47,57,3 5,53,19,51,61,37,33,7,5,11,39,63,27,17,15,23,29,3,21,13,31,25,9,49,33, 
19,29,11,19,27,15,25, 

1,51,61,43,51,59,47,57,3 5,53,19,51,61,37,33,7,5,11,39,63,27,17,15,23,29,3,21,13,31,25,9,4 9,33, 
19,29,11,19,27,15,25, /*end 6 */ 

1,85,67,49,125,25,10 9,43,89,69,113,47,55,97,3,37,83,103,27,13,33,115,41,79,17,2 9,119,75,73,105 
, 7, 5 9, 65, 21, 3, 113, 61, 89, 45, 107, 

1,85,67,49,125,25,109,43,89,69,113,47,55,97,3,37,83,103,27,13,33,115,41,79,17,2 9,119,75,73,105 
, 7, 5 9, 65, 21, 3, 113, 61, 89, 45, 107, 

1,85,67,4 9,12 5,25,109,43,89,69,113,47,55,97,3,37,83,103,2 7,13,33,115,41,79,17,2 9,119,75,73,105 
, 7, 59, 65, 21, 3, 113, 61, 89, 45, 107, 

1,85,67,4 9,12 5,25,109,43,89,69,113,47,55,97,3,37,83,103,2 7,13,3 3,115,41,7 9,17,2 9,119,75,73,105 
, 7, 59, 65,21, 3 , 113 , 61 , 89 , 45 , 107 , 

1,85,67,49,125,2 5,109,43,89,69,113,47,55,97,3,37,83,103,27,13,33,115,41,79,17,29,119,75,73,105 
, 7, 59, 65, 21, 3, 113, 61, 8 9, 45, 107, 

1,85,67,49,125,25,109,43,89,69,113,47,55,97,3,37,83,103,27,13,33,115,41,79,17,29,119,75,73,105 
, 7, 59, 65, 21, 3, 113, 61, 8 9, 45, 107, 



/*end 4 */ 
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1,85,67,49,125,25,109,43,89,69,113,47,55,97,3,37,83,103,27,13,33,115,41,79,17,29,119,75,73,105 
,7,59, 65,21, 3, 113, 61, 89,45, 107, 

1,85,67,49,125,25,109,43,89,69,113,47,55,97,3,37,83,103,27,13,33,115,41,79,17,29,119,75,73,105 
,7,59, 65,21,3/113, 61, 89,45, 107, 

1,85,67,49,125,25,109,43,89,69,113,47,55,97,3,37,83,103,27,13,33,115,41,79,17,29,119,75,73,105 
,7,59,65,21,3,113,61,89,45,107, /*end 7 */ 

1,255,79,147,141,89,173,43,9,25,115,97,19,97,197,101,255,29,203,65,195,177,105,17,47,169,109,1 
49, 15,213, 135,253, 191, 155, 175, 63,47, 7, 23,39, 

1,255,79,147,141,89,173,43,9,25,115,97,19,97,197,101,255,29,203,65,195,177,105,17,47,169,109,1 
49, 15,213, 135,253, 191, 155, 175,63,47, 7,23,39, 

1,255, 79, 147, 141, 89, 173,43, 9,25, 115, 97, 19, 97, 197, 101,255,2 9,2 03,65, 195, 177, 105, 17,47, 169, 109, 1 
4 9, 15, 213 , 13 5, 253 , 191 , 155 , 175 , 63 , 47 , 7 , 2 3 , 3 9 , 

1,255,79,147,141,89,173,43,9,25,115,97,19,97,197,101,255,29,2 03,65,195,177,105,17,47,169,109,1 
49,15,213,135,253,191,155,175,63,47,7,23,39, 

1,255,79,147,141,89,173,43,9,25,115,97,19,97,197,101,255,29,2 03,65,195,177,105,17,47,169,109,1 
4 9, 15,213, 135,253, 191, 155, 175,63,47,7, 23,39, 

1,255,79,147,141,89,173,43,9,25,115,97,19,97,197,101,255,29,2 03,65,195,177,105,17,47,169,109,1 
4 9, 15,213, 135,253, 191, 155, 175,63,47,7, 23,39, 

1,255,79,147,141,89,173,43,9,25,115,97,19,97,197,101,255,29,2 03,65,195,177,105,17,47,169,109,1 
4 9, 15, 213, 13 5, 253 , 191, 155, 175, 63 , 47, 7, 23 , 39, 

i, 255, 79, 147, 141, 89, 173, 43, 9, 25, 115, 97, 19, 97, 197, 101, 255, 29, 2 03, 65, 195, 177, 105, 17, 47, 169, 109,1 
^9, 15, 213, 135, 253 , 191 , 155 , 175 , 63 ,47, 7, 23, 39, 

i§, 255, 79, 147, 141, 89, 173, 43, 9, 25, 115, 97, 19, 97, 197, 101, 255, 29, 2 03, 65, 195, 177, 105, 17,47,169,109,1 
39,15,213,135,253,191,155,175,63,47,7,23,39, /*end 8 */ 

FT: 

II, 2 57,4 65,43 9,177,321,181,225,23 5,103,411,23 3,5 9,353,32 9,463,3 85,111,475,451,263,19,24 9,2 75,3 6 
M/ 3 93, 167, 333,473,469, 101,21,451,229,247, 297,403,497, 61, 361, 

% 257,465,439,177,321,181,225,235,103,411,23 3,59,353,329,463,385,111,475,451,263,19,249,275,36 
Lg, 3 93, 167, 3 3 3,4 73,4 69, 101,21,451,22 9,247, 2 97,403,4 97, 61, 3 61, 

1^,257,465,43 9,177,321,181,225,235,103,411,23 3,59,353,329,463,3 85,111,475,451,263,19,249,275,36 

3 93, 167, 3 3 3,4 73,4 69, 101,21,451,22 9,2 4 7, 2 97,403,4 97, 61, 3 61, 
Zl, 257,465,43 9,177,321,181,225,235,103,411,23 3,59,353,329,463,385,111,475,451,263,19,249,275,36 
^,393, 167,333,473,4 69,101,21,451,229,247,2 97,403,497,61,361, 

III, 257,465,43 9,177,321,181,225,235,103,411,23 3,59,353,32 9,463,3 85, 111,475,451,2 63,19,249,275,36 
gjj, 3 93, 167, 333, 473, 469, 10 1,21, 451, 229, 247, 2 97, 403, 4 97, 61, 3 61, 

; £, 257,465,439,177,321,181,225,235,103,411,233,59,353,329,463,385,111,475,451,263,19,249,275,36 
% 3 93, 167,333,473,469,101,21,451, 229,247,297,403,497,61,361, 

#,257,4 65,43 9,177,321,181,225,235,103,411,233,59,353,32 9,4 63, 385,111,475,451,263,19,249,275,36 
Ci, 393, 167, 333,473,469, 101,21,451,229,247,297,403,497,61,361, 

1,257,4 65,43 9,177,321,181,225,235,103,411,233,59,353,32 9,4 63,3 85,111,475,451,2 63,19,249,275,36 
9,3 93,167,333,473,469, 101,21,451,22 9,247,297,4 03,4 97,61,3 61, 

1,257,4 65,43 9, 177,321,181,225,235,103,411,233,59,353,329,463,385,111,475,451,263,19,249,275,36 
9,393, 167,333,473,469, 101,21,451,229,247,2 97,403,497,61,361, 
/*end 9 */ 

1,771,721,1013,75 9,83 5,94 9,113,92 9,615,157,3 9,761,169,983,657,647,581,505,83 3,13 9,147,2 03,81,3 
37, 829, 989, 375, 365, 131, 215, 733,451,447, 177, 57,471, 979, 197, 2 51, 

1, 771, 721, 1013, 75 9, 835, 94 9, 113 , 92 9 , 615 , 157 , 3 9 , 761 , 169 , 983 , 657 , 647 , 581 , 505 , 83 3 , 13 9, 147, 203, 81, 3 
37, 82 9, 989, 3 75, 3 65, 131, 215, 73 3, 4 51, 447, 177, 57, 4 71, 97 9, 197, 251, 

1, 771, 721, 1013, 75 9, 835 , 949 , 113 , 92 9 , 615 , 157 , 39 , 761, 169 , 983 , 657 , 647 , 581 , 505 , 833 , 139 , 147 , 2 03 , 81, 3 
37,829,989,375,365, 131,215,733,451,447,177,57,471,979,197,251, 

1, 771, 721, 10 13 ,759, 835, 949, 113 , 92 9 , 615 , 157 , 3 9 , 761 , 169 , 983 , 657 , 647 , 5 8 1 , 505 , 83 3,13 9,147,203,81,3 
3 7,82 9,989,375,3 65,131,215,733,451,447,177,57,471,979,197,251, 

1,771,721,1013,759, 835,949,113,929,615,157,39,761,169,983,657,647,581,505,833,139,147,203,81,3 
3 7,82 9,989,3 75,3 65,131,215,73 3,451,447,177,57,471,979,197,2 51, 

1,771,721,1013,759,83 5,94 9,113,92 9,615,157,3 9,761,169,983,657,647,581,505,833,13 9,147,2 03,81,3 
3 7, 82 9, 98 9, 3 75, 3 65, 13 1 , 2 15 , 733 , 451 , 44 7 , 177 , 57 , 471 , 979 , 197 , 25 1 , 

1,771,721,1013,759,835,949,113,929,615,157,39,761,169,983,657,647,581,505,833,13 9,147,203,81,3 
3 7, 82 9, 98 9, 375, 3 65, 131, 215, 733,451,44 7, 177, 57, 471, 97 9, 197, 251, 
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1,771,721,1013,759, 835,949,113,929,615,157,39,761,169,983,657,647,581,505, 833,139,147,203,81,3 
37, 829, 989,375, 365, 131,215,733,451,447,177,57,471,979, 197,2 51, 

1,771,721,1013,75 9,83 5,949,113,929,615,157,3 9,761,169,983,657,647,581,505,83 3,139,147,203,81,3 
37, 82 9, 989, 3 75; 3 65, 131,215,73 3,451,447, 177,57,471, 979, 197,251, 
/*end 10 */ 

1, 1285, 823, 727, 267, 833,471, 1601, 1341, 913, 1725, 2021, 1905,375, 893, 1599,415, 605, 819, 975, 915, 1715, 
1223, 1367, 663,629, 525,469, 981, 1667,1587, 1251,451,481,721,483, 1209, 1457,415, 1435, 

1, 1285, 823, 727,267, 833,471, 1601, 1341, 913, 1725, 2021,1905,375, 893, 1599,415,605, 819, 975, 915, 1715, 
1223, 13 67, 663, 62 9, 525,4 69, 981, 1667, 1587, 12 51,451,481, 721,483, 1209, 1457,415, 1435, 

1, 12 85, 823, 72 7, 267, 833, 471, 1601, 1341, 913 , 1725 , 2 02 1 , 1905 , 3 75 , 8 93 , 15 99 , 415 , 605 , 819 , 975 , 915 , 1715 , 
1223, 1367,663, 629, 525,469, 981, 1667, 1587, 1251,451,481,721,483, 12 09, 1457,415, 1435, 

1,1285,823,727,267,83 3,471,1601,1341,913,1725,2021,1905,3 75,8 93,1599,415,605,819,975,915,1715, 
1223, 13 67,663, 62 9, 52 5,469, 981, 1667, 1587, 1251,451,4 81,721,4 83, 12 09, 1457,415, 143 5, 

1, 12 85, 823,727,267, 83 3,471, 1601, 1341, 913, 1725,2 021, 1905,3 75, 893,1599,415, 605, 819, 975, 915, 1715, 
1223,1367,663,62 9,525,469,981,1667,1587,1251,451,481,721,483,12 09,1457,415,143 5, 

1, 1285, 823, 727, 267, 833,471, 1601, 1341, 913, 1725, 2021, 1905,375, 893, 1599,415, 605, 819, 975, 915, 1715, 
1223,13 67,663,62 9,525,4 69,981,1667,1587,1251,451,4 81,721,4 83,12 09,1457,415,1435, 

1, 1285, 823,727,267, 833,471, 1601, 1341,913, 1725, 2021, 1905,3 75, 893, 1599,415,605, 819, 975, 915, 1715, 
1223, 13 67,663, 62 9, 525,469,981, 1667, 1587, 1251,451,481,721,483, 1209, 1457,415, 1435, 

1,12 85,823,727,267,833,471,1601,1341,913,1725,2021,1905,375,893,1599,415,605,819,975,915,1715, 
J.223, 1367,663, 62 9,525,469, 981, 1667, 1587, 1251,451,481, 721,483, 12 09, 1457,415, 1435, 

3, 12 85, 823,727,267, 83 3,471, 1601, 1341, 913, 1725,2021,1905,375, 8 93, 1599,415,605, 819, 975, 915, 1715, 
^1122 3,1367,663,62 9,52 5,4 69,981,1667,1587,1251,451,481,721,4 83,12 09,1457,415,1435, 
gfend 11 */ 

F? = 

If, 3 855,4 091, 987,183 9,4033,2515,579,3 863,977,3463,2909,3379,1349,373 9,347,3 87,2 881,2821,1873,19 
^,1929,2389,3251,114 9,243,3609,1131,1701,143,1339,3497,2499,1571,983,4021,1625,3217,1163,2977 

UJ, 3855,4091,987,183 9,4033,2515,579,3 863,977,3463,2909,3379,1349,3739,347,387,2 881,2821,1873,19 
J|p, 1929,2389, 3251, 114 9,243, 3609, 1131, 1701, 143 , 1339 , 3497 , 2499 , 1571 , 983 , 4021 , 1625 , 3217 , 1163 , 2977 

% 3855, 4091, 987, 1839,4033,2515,579,3863, 977,3463,2909,3379, 1349,3739,347,387,2881,2821, 1873,19 
619,1929,2389,3251,114 9,243,3609,1131,1701,143,1339,3497,2499,1571,983,4021,162 5,3217,1163,2977 

m 3855, 4091, 987, 183 9, 4033, 2515, 579 ,3 863, 977, 3463, 2909, 3379, 1349, 3739, 347, 387, 2 881, 2821, 1873, 19 
1929, 23 89, 3251, 1149, 243, 3609, 1131, 1701, 143 , 13 3 9 , 34 97 , 2499 , 1571 , 983 , 4021 , 162 5 , 3217 , 1163 , 2977 

®, 3855, 4091, 987, 183 9,4033,2515,579,3 863,977,3463,2909,337 9,1349,373 9,347,3 87,2 881,2821,1873,19 
[£9, 1929, 23 89, 3251, 114 9, 243, 3609, 1131, 1701, 143, 1339, 3497,2499, 1571 , 983 , 402 1 , 1625 , 3217 , 1163 , 2977 

1,3 855,4 091,987, 1839,4033,2515,579,3863,977,3463,2909,3379,1349,3739,347,387,2881,2821,1873,19 
59, 1929, 2389, 3251, 114 9,243, 3609, 1131, 1701, 143 , 133 9 , 34 97 , 24 99 , 1571 , 983 , 4021 , 1625 , 3217 , 1163 , 2977 

1,3855,4091,987,183 9,4033,2515,579,3 863,977,3463,2909,337 9,1349,3739,347,387,2 881,2821,1873,19 
59,192 9,23 89,3251,1149,243,3 60 9,1131,1701,143,13 3 9,34 97,2499,1571,983,4 021,1625,3217,1163,2 977 

1,3 855,4091,987,183 9,4033,2515,579,3 863,977,3463,2909,3379,1349,3739,347,3 87,2 881,2821,1873,19 
5 9,192 9,2 3 89,3251,114 9,243,3609,1131,1701,143,133 9,3497,24 99,1571,983,4021,1625,3217,1163,2 977 

1,3 855,4 091,987,183 9,4 03 3,2515,57 9,3 863,977,3463,2 909,3 3 79,1349,373 9,347,3 87,2 881,2 821,1873,19 
59,192 9,23 8 9,3251,114 9,243,3 60 9,1131,1701,143,13 3 9,34 97,2499,1571,983,4 021,1625,3217,1163,2 977 
/*end 12 */ 

}; 

if (n < 0) { 

for ( j =12 ; j <=MAXBIT; j ++) iv [1+ j *MAXDIM] =1 ; /* Initialize all direction 

numbers for the first 

coordinate to 1 */ 
for ( j =1 , k=0 ; j <=MAXBIT; j ++ , k+-MAXDIM) iu[j] = &iv [k] ; 
for (k=l;k<=MAXDIM;k++) { 

for (j=l; j<=mdeg [k] ; j++) iu[j][k] <<= (MAXBIT-j); 
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for (j=mdeg [k] +1; j<=MAXBIT; j++) { 
ipp=ip [k] ; 
i = iu[j-mdeg[k] ] [k] ; 
i A = (i >> mdeg [k] ) ; 
for (l=mdeg[k] -1/1>=1/1--) { 
if (ipp & l) i A = iu[j-l] [k] ; 

ipp >>= 1; 

} 

iu[j] [k]=i; 



fac=1.0/(lL << MAXBIT); 
in=0; 

} 

else 

■ { 

/* Check if the (n-1) -th number was generated in the previous call 
to sobol . If not, update in and ix */ 

if(in!=n-l) { 
unsigned long gray; 

□ /* Set ix to 0 */ 

U3 f or (k=l;k< = IMIN(d,MAXDIM) ;k+ + ) ix[k]=0; 
ri in=n- 1 ; 

^ii gray=in A (in>>l) ; /* Find gray code of in */ 
yjfor ( j =1 / j <=MAXBIT; j ++) { 

^ if (gray&l) { /* Only digits which are 1 are used */ 

S] im= (j -1) *MAXDIM; 

hi for (k=l;k<=IMIN(d,MAXDIM) ;k++) ix [k] A = iv[im+k] ; 

2 } 

gray>>=l; 

L } 

u } 

|1| im=in; /* Calculate the next vector in the sequence */ 

ps for ( j =1 ; j <=MAXBIT; j ++) { /* Find the rightmost zero bit */ 

*5 if ( ! (im & 1) ) break; 

im >>= 1; 

yj } 

£0 if (j > MAXBIT) nrerror ("MAXBIT too small in sobseq" ) ; 

im= ( j -1) *MAXDIM; 

for (k=l;k<=IMIN(d,MAXDIM) ;k++) { 
ix[k] A = iv[im+k]; 
x [k-1] =ix [k] *fac; 

} 

in++; 



#undef MAXBIT 
#undef MAXDIM 



• 
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/********************************************************************** 

void halton(int n) 

This is the function halton for generating Halton points. 

It returns the n-th d-dimensional Halton point. The point is implicitly 

returned through the array x. The last two digits of n-1 in base p[j] 

are kept in ql[j] and q2[j] . When both digits become p[j]-l, the 

radical inverse function is computed again. That way the accumulation 

of round-off error is avoided. In practice, there are not any upper bounds 

on the values of d and n. 



********************* 



*************************************************/ 



extern int d; /* actual dimension of the points */ 

extern int *ql / *q2; /* ql[j] is the last digit of n-1 in base p[j], q2[j] is 

the digit before the last one */ 

extern double *x; /* This will contain the Halton point */ 

extern int *p; /* the first d prime numbers */ 

extern int *p_l; /* first d prime numbers minus 1 */ 

extern double *incrl, *incr2 ; /* incrl[j] is l/p[j] and 1/ (p [ j ] *p [ j ] ) */ 

Rouble find_fi(int p, int n) ; /* See below */ 

yyi 

[jfbid halton (int n) 

^ double a; 
«jj int j ,nn; 

~™4 static int ins_n; /* The default value of ins__n is 0 */ 

hi /* Check if the (n-1) -th number was generated in the previous call 
^" to halton. If not, update ql, q2 , and x */ 

D if (ins n ! =n-l) 

ru { " 

ins_n=n-l; 
% for (j=0; j<d; j++) 

%y ql [ j ] =ins_n%p [ j ] ; 
08 q2 [ j ] = (ins_n/p [ j ] ) %p [ j ] ; 
x [ j ] =f ind_f i (p [j ] , ins_n) ; 

} 

ins n++; 

} 

for (j=0; j<d; j++) 
{ 

if (ql[j]<p_l[jl) 

{ 

/* It is easy to update when the last digit is less than p[j]-l */ 
ql [j] ++; 

x[j]=x[j]+incrl[j] ; 

} 

else if (q2[j]<p_l[j]) 

{ 

/* This is the case when the last digit is ptj]-l and the digit 
before the last one is less than p[j]-l */ 



ql[j]=0; 
q2 [j] ++; 



■25- 



# 



# 
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} 

} 



x t j ] =x [ j ] +incrl [ j ] +incr2 [ j ] -1.0; 
else 



/* This is the case when the last digit is p[j]-l and the digit 
before the last one is also p[j]-l */ 

ql[j]=0; 
q2 [j]=0; 

nn=n/ (p [j ] *p [j ] ) ; 
a=nn%p [j] / 
nn=nn/p [j ] ; 

if (nn) x [ j ] = (a+f ind_f i (p [ j ] , nn) ) *incr2 [ j ] *incrl [ j ] / 
else x [ j ] =a*incr2 [ j ] *incrl [ j ] ; 



} 



find_fi(int p, int n) 

*rfThis returns the radical inverse function fi(p,n) at n for the prime p. 
*********************^ 

i&uble find fi(int p, int n) 

If. 

' T -j int p2 , nn; 
LU double s,fi,incr; 

K incr=1.0/p; 
_L ; p2=p*p; 
^ nn=n/p2 ; 

i 

Qfi^O.O; 

s = incr ; 
^ while (nn > 0) 

6B fi+= (nn%p)*s; 
nn-nn/p; 
s*=incr ; 

} 

/* The two largest components of fi are added later to avoid possible 

loss of precision */ 
fi=((n/p%p)+ fi)/ p2; 
fi+=(n%p)/ (double) p; 
return fi; 

} 



